Prep pCVP2-BRI1-mCitrine seu object with metadata¶

In [1]:
library(tidyverse)
library(Seurat)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
library(muscat)
library(purrr)
library(limma)
library(scran)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.1     ✔ tibble    3.1.8
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching SeuratObject


Attaching package: ‘cowplot’


The following object is masked from ‘package:lubridate’:

    stamp


Loading required package: grid

========================================
ComplexHeatmap version 2.14.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite either one:
- Gu, Z. Complex Heatmap Visualization. iMeta 2022.
- Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
    genomic data. Bioinformatics 2016.


The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================


========================================
circlize version 0.4.15
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
  in R. Bioinformatics 2014.

This message can be suppressed by:
  suppressPackageStartupMessages(library(circlize))
========================================


Loading required package: SingleCellExperiment

Loading required package: SummarizedExperiment

Loading required package: MatrixGenerics

Loading required package: matrixStats


Attaching package: ‘matrixStats’


The following object is masked from ‘package:dplyr’:

    count



Attaching package: ‘MatrixGenerics’


The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars


Loading required package: GenomicRanges

Loading required package: stats4

Loading required package: BiocGenerics


Attaching package: ‘BiocGenerics’


The following object is masked from ‘package:limma’:

    plotMA


The following objects are masked from ‘package:lubridate’:

    intersect, setdiff, union


The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union


The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs


The following objects are masked from ‘package:base’:

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
    Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min


Loading required package: S4Vectors


Attaching package: ‘S4Vectors’


The following objects are masked from ‘package:lubridate’:

    second, second<-


The following objects are masked from ‘package:dplyr’:

    first, rename


The following object is masked from ‘package:tidyr’:

    expand


The following objects are masked from ‘package:base’:

    expand.grid, I, unname


Loading required package: IRanges


Attaching package: ‘IRanges’


The following object is masked from ‘package:lubridate’:

    %within%


The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice


The following object is masked from ‘package:purrr’:

    reduce


Loading required package: GenomeInfoDb

Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.



Attaching package: ‘Biobase’


The following object is masked from ‘package:MatrixGenerics’:

    rowMedians


The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians



Attaching package: ‘SummarizedExperiment’


The following object is masked from ‘package:SeuratObject’:

    Assays


The following object is masked from ‘package:Seurat’:

    Assays


Loading required package: scuttle

In [2]:
library(future)
#for 200gb ram 
options(future.globals.maxSize = 200000 * 1024^2)
In [3]:
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Stream 8

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/tmn23/miniconda3/envs/muscat/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    grid      stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] future_1.31.0               scran_1.26.0               
 [3] scuttle_1.8.0               SingleCellExperiment_1.20.0
 [5] SummarizedExperiment_1.28.0 Biobase_2.58.0             
 [7] GenomicRanges_1.50.0        GenomeInfoDb_1.34.8        
 [9] IRanges_2.32.0              S4Vectors_0.36.0           
[11] BiocGenerics_0.44.0         MatrixGenerics_1.10.0      
[13] matrixStats_0.63.0          limma_3.54.0               
[15] muscat_1.12.0               ggrepel_0.9.3              
[17] gprofiler2_0.2.1            GeneOverlap_1.34.0         
[19] circlize_0.4.15             ComplexHeatmap_2.14.0      
[21] cowplot_1.1.1               SeuratObject_4.1.3         
[23] Seurat_4.3.0                lubridate_1.9.2            
[25] forcats_1.0.0               stringr_1.5.0              
[27] dplyr_1.1.0                 purrr_1.0.1                
[29] readr_2.1.4                 tidyr_1.3.0                
[31] tibble_3.1.8                ggplot2_3.4.1              
[33] tidyverse_2.0.0            

loaded via a namespace (and not attached):
  [1] pbdZMQ_0.3-9              scattermore_0.8          
  [3] bit64_4.0.5               irlba_2.3.5.1            
  [5] DelayedArray_0.24.0       data.table_1.14.8        
  [7] KEGGREST_1.38.0           RCurl_1.98-1.10          
  [9] doParallel_1.0.17         generics_0.1.3           
 [11] ScaledMatrix_1.6.0        RhpcBLASctl_0.23-42      
 [13] RSQLite_2.2.20            RANN_2.6.1               
 [15] bit_4.0.5                 tzdb_0.3.0               
 [17] spatstat.data_3.0-0       httpuv_1.6.9             
 [19] viridis_0.6.2             hms_1.1.2                
 [21] evaluate_0.20             promises_1.2.0.1         
 [23] fansi_1.0.4               progress_1.2.2           
 [25] caTools_1.18.2            igraph_1.3.5             
 [27] DBI_1.1.3                 geneplotter_1.76.0       
 [29] htmlwidgets_1.6.1         spatstat.geom_3.0-6      
 [31] ellipsis_0.3.2            backports_1.4.1          
 [33] annotate_1.76.0           aod_1.3.2                
 [35] deldir_1.0-6              sparseMatrixStats_1.10.0 
 [37] vctrs_0.5.2               ROCR_1.0-11              
 [39] abind_1.4-5               cachem_1.0.6             
 [41] withr_2.5.0               progressr_0.13.0         
 [43] sctransform_0.3.5         prettyunits_1.1.1        
 [45] goftest_1.2-3             cluster_2.1.4            
 [47] IRdisplay_1.1             lazyeval_0.2.2           
 [49] crayon_1.5.2              genefilter_1.80.0        
 [51] spatstat.explore_3.0-6    edgeR_3.40.0             
 [53] pkgconfig_2.0.3           nlme_3.1-162             
 [55] vipor_0.4.5               blme_1.0-5               
 [57] rlang_1.0.6               globals_0.16.2           
 [59] lifecycle_1.0.3           miniUI_0.1.1.1           
 [61] rsvd_1.0.5                polyclip_1.10-4          
 [63] lmtest_0.9-40             Matrix_1.5-3             
 [65] IRkernel_1.3.2            boot_1.3-28.1            
 [67] zoo_1.8-11                base64enc_0.1-3          
 [69] beeswarm_0.4.0            ggridges_0.5.4           
 [71] GlobalOptions_0.1.2       png_0.1-8                
 [73] viridisLite_0.4.1         rjson_0.2.21             
 [75] bitops_1.0-7              KernSmooth_2.23-20       
 [77] Biostrings_2.66.0         blob_1.2.3               
 [79] DelayedMatrixStats_1.20.0 shape_1.4.6              
 [81] parallelly_1.34.0         spatstat.random_3.1-3    
 [83] beachmat_2.14.0           scales_1.2.1             
 [85] memoise_2.0.1             magrittr_2.0.3           
 [87] plyr_1.8.8                ica_1.0-3                
 [89] gplots_3.1.3              zlibbioc_1.44.0          
 [91] compiler_4.2.2            dqrng_0.3.0              
 [93] RColorBrewer_1.1-3        clue_0.3-64              
 [95] lme4_1.1-31               DESeq2_1.38.0            
 [97] fitdistrplus_1.1-8        cli_3.6.0                
 [99] XVector_0.38.0            lmerTest_3.1-3           
[101] listenv_0.9.0             patchwork_1.1.2          
[103] pbapply_1.7-0             TMB_1.9.2                
[105] MASS_7.3-58.2             tidyselect_1.2.0         
[107] stringi_1.7.12            BiocSingular_1.14.0      
[109] locfit_1.5-9.7            tools_4.2.2              
[111] timechange_0.2.0          future.apply_1.10.0      
[113] parallel_4.2.2            uuid_1.1-0               
[115] bluster_1.8.0             foreach_1.5.2            
[117] metapod_1.6.0             gridExtra_2.3            
[119] Rtsne_0.16                digest_0.6.31            
[121] shiny_1.7.4               Rcpp_1.0.10              
[123] broom_1.0.3               later_1.3.0              
[125] RcppAnnoy_0.0.20          httr_1.4.4               
[127] AnnotationDbi_1.60.0      Rdpack_2.4               
[129] colorspace_2.1-0          XML_3.99-0.13            
[131] tensor_1.5                reticulate_1.28          
[133] splines_4.2.2             statmod_1.5.0            
[135] uwot_0.1.14               spatstat.utils_3.0-1     
[137] scater_1.26.0             sp_1.6-0                 
[139] plotly_4.10.1             xtable_1.8-4             
[141] jsonlite_1.8.4            nloptr_2.0.3             
[143] R6_2.5.1                  pillar_1.8.1             
[145] htmltools_0.5.4           mime_0.12                
[147] glue_1.6.2                fastmap_1.1.0            
[149] minqa_1.2.5               BiocParallel_1.32.5      
[151] BiocNeighbors_1.16.0      codetools_0.2-19         
[153] utf8_1.2.3                lattice_0.20-45          
[155] spatstat.sparse_3.0-0     numDeriv_2016.8-1.1      
[157] pbkrtest_0.5.2            ggbeeswarm_0.7.1         
[159] leiden_0.4.3              gtools_3.9.4             
[161] survival_3.5-3            glmmTMB_1.1.5            
[163] repr_1.1.6                munsell_0.5.0            
[165] GetoptLong_1.0.5          GenomeInfoDbData_1.2.9   
[167] iterators_1.0.14          variancePartition_1.28.0 
[169] reshape2_1.4.4            gtable_0.3.1             
[171] rbibutils_2.2.13         
In [4]:
rc.integrated <- readRDS("../../CheWei/scRNA-seq/Integrated_Objects/rc.integrated_11S_CVP_BRI1_seu3_20230315.rds")
In [5]:
rc.integrated
An object of class Seurat 
71788 features across 72381 samples within 3 assays 
Active assay: integrated (17520 features, 17520 variable features)
 2 other assays present: RNA, SCT
 4 dimensional reductions calculated: pca, umap, umap_3D, umap_2D
In [6]:
# remove sc_12 dc1 dc2
rc.integrated <- subset(rc.integrated, subset = orig.ident %in% c("briT", "briTR","sc_130", "sc_131", "sc_132", "sc_134", "sc_135", "sc_136"))
In [7]:
table(rc.integrated$orig.ident)
  briT  briTR sc_130 sc_131 sc_132 sc_134 sc_135 sc_136 
  7483   8318   6589   7615   6550   7745   5038   6089 
In [8]:
rc.integrated$geno <- rc.integrated$orig.ident
rc.integrated$geno <- gsub("briT","bri1_T",rc.integrated$geno) 
rc.integrated$geno <- gsub("bri1_TR","pCVP2_BRI1_Citrine_bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_130","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_131","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_132","pCVP2_BRI1_Citrine_bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_134","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_135","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_136","pCVP2_BRI1_Citrine_bri1_T",rc.integrated$geno)
In [9]:
rc.integrated$source <- rc.integrated$orig.ident
rc.integrated$source <- gsub("briT","Graeff-2021",rc.integrated$source) 
rc.integrated$source <- gsub("bri1_TR","Graeff-2021",rc.integrated$source)
rc.integrated$source <- gsub("Graeff-2021R","Graeff-2021",rc.integrated$source)
rc.integrated$source <- gsub("sc_130","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_131","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_132","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_134","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_135","Benfey",rc.integrated$source)
rc.integrated$source <- gsub("sc_136","Benfey",rc.integrated$source)
In [10]:
rc.integrated$rep <- rc.integrated$orig.ident
rc.integrated$rep <- gsub("briT","1",rc.integrated$rep)
rc.integrated$rep <- gsub("bri1_TR","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_130","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_131","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_132","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_134","3",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_135","3",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_136","3",rc.integrated$rep)
In [11]:
table(rc.integrated$orig.ident, rc.integrated$geno)
        
         bri1_T pCVP2_BRI1_Citrine_bri1_T   WT
  briT     7483                         0    0
  briTR       0                      8318    0
  sc_130      0                         0 6589
  sc_131   7615                         0    0
  sc_132      0                      6550    0
  sc_134      0                         0 7745
  sc_135   5038                         0    0
  sc_136      0                      6089    0
In [30]:
table(rc.integrated$orig.ident, rc.integrated$source)
        
         Benfey Graeff-2021
  briT        0        7483
  briTR       0        8318
  sc_130   6589           0
  sc_131   7615           0
  sc_132   6550           0
  sc_134   7745           0
  sc_135   5038           0
  sc_136   6089           0
In [29]:

In [14]:
feature_names <- read_tsv("../data/features.tsv.gz", col_names = c("AGI", "Name", "Type")) %>%
  select(-Type) %>%
  distinct()
Rows: 32833 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): AGI, Name, Type

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
In [15]:
rc.integrated$geno <- factor(rc.integrated$geno, levels=c("WT", "bri1_T", "pCVP2_BRI1_Citrine_bri1_T"))

Cell and developmental stage metadata¶

  • Developmental stage: time_zone
  • Cell type:cell_type
  • Combination of cell type and developmental stage: time_zone_cell_type
  • Combination of cell type and developmental stage with cell subtypes (not used): time_zone_cell_subtypes
In [16]:
rc.integrated$cell_type <- rc.integrated$celltype.anno.Li.crude
rc.integrated$time_zone <- rc.integrated$time.anno
rc.integrated$time_zone_cell_type <- rc.integrated$time.celltype.anno.Li.crude
In [17]:
table(rc.integrated$orig.ident, rc.integrated$cell_type)
        
         Atrichoblast Columella Cortex Endodermis Lateral Root Cap Pericycle
  briT            640       403    332        401             2368      1750
  briTR           700       637    228        442             4148      1167
  sc_130         1445       895    972        725             1058       143
  sc_131         1441       672    777        490             1800       812
  sc_132         1089       474    720        544             1597       554
  sc_134         1171       582    906        573             2101       633
  sc_135         1055       294    544        480              991       428
  sc_136         1176       649    878        617             1240       247
        
         Phloem Procambium Quiescent Center Trichoblast Xylem
  briT      204        815                7         200   363
  briTR     227        358               13         148   250
  sc_130     52        127                0        1000   172
  sc_131    130        320                0         900   273
  sc_132    138        288                0         916   230
  sc_134    123        348                0        1041   267
  sc_135     88        260                0         707   191
  sc_136     78        186                0         871   147
In [18]:
table(rc.integrated$time.anno)
         Distal Columella   Distal Lateral Root Cap                Elongation 
                     3943                      5190                     15349 
               Maturation                  Meristem        Proximal Columella 
                    10718                      9559                       633 
Proximal Lateral Root Cap 
                    10035 
In [19]:
# Plot celltype annotation Li
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400D3", "#DCD0FF","#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#E6194B", "#DD77EC", "#9A6324", "#FFE119", "#FF9900", "#FFD4E3", "#9A6324", "#DDAA6F", "#EEEEEE")
rc.integrated$cell_type <- factor(rc.integrated$cell_type, levels = order[sort(match(unique(rc.integrated$cell_type),order))])
color <- palette[sort(match(unique(rc.integrated$cell_type),order))]


options(repr.plot.width=15, repr.plot.height=7)

(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'geno', pt.size = 0.75, ncol=3))
In [27]:
options(repr.plot.width=15, repr.plot.height=7)
FeaturePlot(rc.integrated, features="BRI1-mCitrine", split.by = "geno", order=T)
Warning message in FeaturePlot(rc.integrated, features = "BRI1-mCitrine", split.by = "geno", :
“All cells have the same value (0) of BRI1-mCitrine.”
In [22]:
options(repr.plot.width=30, repr.plot.height=7)

(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'orig.ident', pt.size = 0.75, ncol=8))
In [23]:
DefaultAssay(rc.integrated) <- "SCT"
FeaturePlot(rc.integrated, features="BRI1-mCitrine", split.by = "orig.ident", order=T)
Warning message in FeaturePlot(rc.integrated, features = "BRI1-mCitrine", split.by = "orig.ident", :
“All cells have the same value (0) of BRI1-mCitrine.”
Warning message in FeaturePlot(rc.integrated, features = "BRI1-mCitrine", split.by = "orig.ident", :
“All cells have the same value (0) of BRI1-mCitrine.”
In [31]:
saveRDS(rc.integrated, file = "/hpc/group/pbenfeylab/CheWei/scRNA-seq/Integrated_Objects/rc.integrated_8S_CVP_BRI1_seu3_annotated_20230316.rds")
In [ ]: